setwd("~/ProfessionAI - Data Science/03_Statistica inferenziale")
neonati = read.csv("neonati.csv", stringsAsFactors = T)Modello Statistico per la Previsione del Peso Neonatale
1 Introduzione
L’analisi nasce dall’esigenza di avere una previsione anticipata del peso del neonato, al fine di una più corretta pianificazione delle risorse da allocare in caso di peso inferiore alle aspettative.
Il peso è difatti una caratteristica indicativa e di grande rilievo per la salute del neonato e in caso sia troppo basso, devono essere predisposte strumentazioni e persone per effettuare attività ben precise, atte a riportare il peso entro i parametri corretti.
Il progetto si inserisce all’interno di un contesto di crescente attenzione verso la prevenzione delle complicazioni neonatali. La possibilità di prevedere il peso alla nascita dei neonati rappresenta un’opportunità fondamentale per migliorare la pianificazione clinica e ridurre i rischi associati a nascite problematiche, come parti prematuri o neonati con basso peso. Di seguito, i principali benefici che questo progetto porterà all’azienda e al settore sanitario:
Miglioramento delle previsioni cliniche:
- Il peso del neonato è un indicatore chiave della sua salute. Avere un modello predittivo accurato consente al personale medico di intervenire tempestivamente in caso di anomalie, riducendo le complicazioni perinatali come le difficoltà respiratorie o l’ipoglicemia.
Ottimizzazione delle risorse ospedaliere:
- Sapere in anticipo quali neonati potrebbero avere bisogno di cure intensive aiuta a organizzare le risorse umane e tecnologiche degli ospedali in modo efficiente. Questo si traduce in una riduzione dei costi operativi e una migliore pianificazione dell’utilizzo delle unità di terapia intensiva neonatale (TIN).
Prevenzione e identificazione dei fattori di rischio:
- Il modello potrà evidenziare i fattori che maggiormente influenzano negativamente il peso del neonato (come il fumo materno, gravidanze multiple o età avanzata della madre). Queste informazioni sono preziose per la prevenzione e la gestione personalizzata delle gravidanze, permettendo interventi proattivi in caso di rischio elevato.
Valutazione delle pratiche ospedaliere:
- Attraverso un’analisi comparativa tra i tre ospedali coinvolti, l’azienda potrà identificare eventuali differenze nei risultati clinici, come una maggiore incidenza di parti cesarei in una determinata struttura. Ciò consente di monitorare la qualità delle pratiche e armonizzare i protocolli tra i diversi centri ospedalieri, migliorando la coerenza delle cure.
Supporto alla pianificazione strategica:
- L’analisi dei dati e le previsioni possono essere utilizzate per prendere decisioni informate non solo a livello clinico ma anche strategico. L’azienda potrà sfruttare queste informazioni per implementare nuove politiche di salute pubblica, garantendo un impatto positivo sui tassi di mortalità e morbilità neonatale.
2 Analisi descrittiva
In questa sezione si utilizzeranno i più noti indicatori statistici per descrivere il dataset utilizzato al fine di:
conoscerne adeguatamente le variabili,
individuare eventuali valori anomali e inquadrarli
comprendere la struttura generale dei dati
2.1 Il dataset
2.2 Le librerie
Di seguito le librerie utilizzate per condurre l’analisi. Sostanzialmente tutto il pacchetto tidyverse per analisi, manipolazione dei dati, nonchè creazioni di grafici.
Le librerie car, MASS e lmtest tutte dedicate ai test sul modello statistico definito.
library(ggplot2)
library(dplyr)
library(patchwork)
library(moments)
library(tidyr)
library(car)
library(MASS)
library(lmtest)2.3 Le variabili
I dati analizzati sono stati raccolti da 3 ospedali differenti, in ciascuno dei quali le variabili registrate sono state 8, con 2500 osservazioni totali.
Queste sono le 8 variabili studiate all’interno del campione
colnames(neonati) [1] "Anni.madre" "N.gravidanze" "Fumatrici" "Gestazione" "Peso"
[6] "Lunghezza" "Cranio" "Tipo.parto" "Ospedale" "Sesso" Di seguito una loro breve descrizione:
Anni.madre: Misura dell’età in anni.
N.gravidanze: Quante gravidanze ha avuto la madre.
Fumatrici: Un indicatore binario (0=non fumatrice, 1=fumatrice).
Gestazione: Numero di settimane di gestazione.
Peso: Peso alla nascita in grammi.
Lunghezza: lunghezza del neonato in mm, misurabile anche durante la gravidanza tramite ecografie.
Cranio: diametro craniale, misurabile anche durante la gravidanza tramite ecografie.
Tipo.parto: Naturale o cesareo.
Ospedale: Ospedale di nascita 1, 2 o 3.
Sesso: Maschio (M) o femmina (F).
Di seguito invece una visione di sintesi dei quantili e della media di ciascuna variabile nel dataset
summary(neonati) Anni.madre N.gravidanze Fumatrici Gestazione
Min. : 0.00 Min. : 0.0000 Min. :0.0000 Min. :25.00
1st Qu.:25.00 1st Qu.: 0.0000 1st Qu.:0.0000 1st Qu.:38.00
Median :28.00 Median : 1.0000 Median :0.0000 Median :39.00
Mean :28.16 Mean : 0.9812 Mean :0.0416 Mean :38.98
3rd Qu.:32.00 3rd Qu.: 1.0000 3rd Qu.:0.0000 3rd Qu.:40.00
Max. :46.00 Max. :12.0000 Max. :1.0000 Max. :43.00
Peso Lunghezza Cranio Tipo.parto Ospedale Sesso
Min. : 830 Min. :310.0 Min. :235 Ces: 728 osp1:816 F:1256
1st Qu.:2990 1st Qu.:480.0 1st Qu.:330 Nat:1772 osp2:849 M:1244
Median :3300 Median :500.0 Median :340 osp3:835
Mean :3284 Mean :494.7 Mean :340
3rd Qu.:3620 3rd Qu.:510.0 3rd Qu.:350
Max. :4930 Max. :565.0 Max. :390 Balza subito all’occhio che ci sia qualcosa che non vada nella variabile Anni.madre, poichè presenta un valore minimo di 0, cosa ovviamente impossibile.
Si effettua quindi una verifica più puntuale per individuare altri valori chiaramente anonimi, filtrando il dataset per la variabile anni.madre < 15
neonati %>%
filter(Anni.madre < 15) %>%
arrange(Anni.madre)| Anni.madre | N.gravidanze | Fumatrici | Gestazione | Peso | Lunghezza | Cranio | Tipo.parto | Ospedale | Sesso |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 0 | 39 | 3060 | 490 | 330 | Nat | osp3 | M |
| 1 | 1 | 0 | 41 | 3250 | 490 | 350 | Nat | osp2 | F |
| 13 | 0 | 0 | 38 | 2760 | 470 | 325 | Nat | osp2 | F |
| 14 | 1 | 0 | 39 | 3510 | 490 | 365 | Nat | osp2 | M |
| 14 | 0 | 0 | 39 | 3550 | 500 | 355 | Ces | osp1 | M |
Si evince che, pur essendo valori di età bassi, 13 e 14 non possono essere scartati perchè comunque verosimili. Lo stesso non si può ovviamente dire per 1 e 0.
Probabilmente si tratta di errori di battitura che però potrebbero influenzare negativamente le conclusioni che si prenderanno. Si decide dunque di eliminare queste 2 osservazioni dal dataset
neonati.clean = neonati %>%
filter(Anni.madre > 1)2.4 Gli indici
Di seguito verranno calcolati gli indici di variabilità e di forma, così da avere un’idea della possibile distribuzione delle variabili e soffermarci su envtuali valori degni di nota.
Resteranno fuori da questa analisi quantitativa, le variabili qualitative come Sesso e Ospedale.
2.4.1 Indici di variabilità
Gli indici di variabilità forniscono un secondo livello di informazioni sulle variabili, in particolare indicano quanto esse siano disperse o concentrate attorno al valore medio.
Di seguito una tabella che raccoglie, per ogni variabile quantitativa, gli indici di variabilità corrispettivi.
attach(neonati.clean)
variance = round(c(var(Anni.madre),
var(N.gravidanze),
var(Gestazione),
var(Peso),
var(Lunghezza),
var(Cranio)),2)
standard_dev = round(sqrt(variance),2)
IQR = round(c(IQR(Anni.madre),
IQR(N.gravidanze),
IQR(Gestazione),
IQR(Peso),
IQR(Lunghezza),
IQR(Cranio)),2)
variation_index = round(c(sd(Anni.madre)/mean(Anni.madre),
sd(N.gravidanze)/mean(N.gravidanze),
sd(Gestazione)/mean(Gestazione),
sd(Peso)/mean(Peso),
sd(Lunghezza)/mean(Lunghezza),
sd(Cranio)/mean(Cranio)),2)
VARIATION_GLOBAL = data.frame(IQR, variance, standard_dev, variation_index)
rows = c("anni_madre", "n_gravidanze", "gestazione", "peso", "lunghezza", "cranio")
row.names(VARIATION_GLOBAL) = rows
VARIATION_GLOBAL| IQR | variance | standard_dev | variation_index | |
|---|---|---|---|---|
| anni_madre | 7 | 27.22 | 5.22 | 0.19 |
| n_gravidanze | 1 | 1.64 | 1.28 | 1.30 |
| gestazione | 2 | 3.49 | 1.87 | 0.05 |
| peso | 630 | 275865.90 | 525.23 | 0.16 |
| lunghezza | 30 | 693.21 | 26.33 | 0.05 |
| cranio | 20 | 269.93 | 16.43 | 0.05 |
2.4.2 Indici di forma
Gli indici di forma forniscono informazioni sulla forma della distribuzione, indicando quanto essa sia asimmetrica, con il vertice della campana dunque spostato verso sinistra o destra rispetto a una gaussiana standard, oppure quanto tale vertice si alzi o abbassi rispetto a una gaussiana standard.
Di seguito una tabella con l’indice di asimmetria e curtosi per ogni variabile
Skewness = round(c(skewness(Anni.madre),
skewness(N.gravidanze),
skewness(Gestazione),
skewness(Peso),
skewness(Lunghezza),
skewness(Cranio)),2)
Kurtosis = round(c(kurtosis(Anni.madre),
kurtosis(N.gravidanze),
kurtosis(Gestazione),
kurtosis(Peso),
kurtosis(Lunghezza),
kurtosis(Cranio)),2)
SHAPE_GLOBAL = data.frame(Skewness,Kurtosis)
row.names(SHAPE_GLOBAL) = rows
SHAPE_GLOBAL| Skewness | Kurtosis | |
|---|---|---|
| anni_madre | 0.15 | 2.89 |
| n_gravidanze | 2.51 | 13.98 |
| gestazione | -2.07 | 11.26 |
| peso | -0.65 | 5.03 |
| lunghezza | -1.51 | 9.48 |
| cranio | -0.79 | 5.94 |
2.4.3 Commenti sugli indici
Possiamo senza dubbio concludere che le variabili osservate non presentano un’eccessiva varibilità. Confrontando l’indice di variazione infatti, vale a dire la media divisa per la deviazione standard, solo il numero di gravidanze sembra mostrare un particolare scostamento, con la deviazione standard il 30% più grande della media.
Il numero di gravidanze, insieme alle settimane di gestazione, presentano invece la maggiore asimmetria. In particolare:
- il numero di gravidanze ha un’asimmetria positiva, verso sinistra dunque. Questo ci dice che, tendenzialmente, il nostro campione è composto da mamme con pochi parti alle spalle. La coda destra è comunque lunga, ad indicare che, anche se in poche occorrenze, c’è chi era alla 13ma gestazione.
- Il numero di settimane di gestazione invece ha asimmetria negativa, vale a dire verso destra. Il risultato non stupisce, poichè il numero tipico di settimane di gestazione è 40. Essendo un valore “massimo”, non può generare una gaussiana standard, altrimenti si avrebbe il 50% di probabilità di avere parti oltre la 40esima. Va posta attenzione sul fatto che però ci sia una coda sinistra pronunciata, con casi di nascita addirittura alla 25ma settimana. Questo potrebbe complicare le cose in fase di modellazione.
La Kurtosi invece ci indica la chiara tendenza di tutte le variabili a essere lepticurtiche, cioè con una “campana” più stretta e alta rispetto a una normale standard. Proprio come per l’asimmetria, le variabili decisamente più lepticurtiche, sono proprio le settimane di gestazione e il numero di gravidanze pregresse.
3 Test statistici sul campione
In questo paragrafo verranno sondate alcune ipotesi che i dati lasciano supporre. Questo passo si rende necessario perchè, nonostante i dati del campione possano puntare in una certa direzione, bisogna sempre ricordare che, anche se numeroso, si tratta sempre di un campione.
Pertanto, prima di generalizzare quello che il campione mostra, si deve cercare di capire se ci sia margine o meno per effettuare inferenza.
3.1 Differenze fisiche tra i sessi
La prima domanda cui vogliamo dare risposta è: esiste differenza nelle misure antropometriche (lunghezza, peso e diametro craniale) tra i due sessi?
Si pone nuovamente l’attenzione sul fatto che la domanda sia generale. Non si sta chiedendo se esistano nel campione, delle differenze tra sessi. Per quello un grafico sarebbe più che sufficiente.
Quello che si sta chiedendo è se si possa estendere a una popolazione quanto visto sul campione.
Innanzitutto, vediamo quale sia l’andamento nel campione, relativamente alle misure antropometriche, per i due sessi.
I 3 grafici mostrano che, nel campione, esiste una differenza tra maschi e femmine, per tutte le misure antropometriche, laddove i maschi le hanno tendenzialmente maggiori rispetto alle femmine.
Degno di nota il fatto che esistono molti valori non racchiusi dal boxplot e che dovranno essere studiati in fase di definizione del modello lineare.
Testiamo ora il seguente sistema di ipotesi:
\[ \begin{cases} H_{0} : \mu_{M} - \mu_{F} = 0, & \mbox{La media per i maschi } \mu_{M}\mbox{ coincide coincide con la media per le femmine } \mu_{F}\\ H_{1} : \mu_{M} - \mu_{F} \neq 0, & \mbox{Le medie tra maschi e femmine sono diverse } \end{cases} \]
Di volta in volta assumeremo la media per i 3 indicatori antropometrici: lunghezza, peso e diametro craniale.
Si userà il test t, essendo il più adatto per confronti di medie tra due gruppi. Il livello di significatività (cioè la probabilità di rifiutare l’ipotesi nulla quando in realtà è vera) sarà fissato al 5%.
3.1.1 Differenza di peso tra i sessi
Il test verificherà che la media di peso dei maschi, sottratta alla media di peso delle femmine sia uguale a 0, vale a dire che le due media siano identiche.
Il test assume che le due medie siano uguali. Se questo fosse vero, la differenza tra le due, estraendo dei campioni, si distribuirebbe secondo una curva t di Student.
t.test(Peso~Sesso, data = neonati.clean)
Welch Two Sample t-test
data: Peso by Sesso
t = -12.115, df = 2488.7, p-value < 2.2e-16
alternative hypothesis: true difference in means between group F and group M is not equal to 0
95 percent confidence interval:
-287.4841 -207.3844
sample estimates:
mean in group F mean in group M
3161.061 3408.496 Il test mostra un p-value praticamente a 0. Questo valore indica la probabilità di ottenere una differenza tra le medie uguale rispetto a quella vista in questo campione o peggio, se fosse vera l’ipotesi nulla.
Questo vuol dire che la probabilità che il nostro campione presenti una differenza tra le medie così diversa da 0 rasenta l’impossibilità. Talmente improbabile che diventa assolutamente più probabile che sia sbagliata l’ipotesi per cui le 2 medie sono coincidenti.
Ergo, il test dimostra che le due medie sono diverse tra loro in maniera statisticamente molto significativa.
Possiamo pertanto estendere a tutta la popolazione l’affermazione che: c’è differenza di peso tra i neonati a dipendere dal sesso.
3.1.2 Differenza di diametro craniale tra i sessi
t.test(Cranio~Sesso, data = neonati.clean)
Welch Two Sample t-test
data: Cranio by Sesso
t = -7.4366, df = 2489.4, p-value = 1.414e-13
alternative hypothesis: true difference in means between group F and group M is not equal to 0
95 percent confidence interval:
-6.110504 -3.560417
sample estimates:
mean in group F mean in group M
337.6231 342.4586 Anche in questo caso abbiamo un p-value vicinissimo allo 0. Si rifiuta pertanto l’ipotesi nulla di uguaglianza tra le due medie e possiamo concludere che esiste differenza, tra maschi e femmine, nel diametro craniale.
3.1.3 Differenza di lunghezza tra i sessi
t.test(Lunghezza~Sesso, data = neonati.clean)
Welch Two Sample t-test
data: Lunghezza by Sesso
t = -9.5823, df = 2457.3, p-value < 2.2e-16
alternative hypothesis: true difference in means between group F and group M is not equal to 0
95 percent confidence interval:
-11.939001 -7.882672
sample estimates:
mean in group F mean in group M
489.7641 499.6750 Anche in questo caso il p-value è vicino allo 0 e si rifiuta l’ipotesi nulla. Anche per la lunghezza possiamo affermare che esiste differenza tra maschi e femmine in generale.
3.2 Peso e lunghezza del campione vs popolazione
Vogliamo verificare che il dato medio di peso e lunghezza del nostro campione, sia allineato ai valori medi della popolazione. Il livello di significatività lo si pone pari al 5%.
\[ \begin{cases} H_{0} : \mbox{La media del campione coincide con la media della popolazione } \\ H_{1} : \mbox{La media del campione è diversa dalla media della popolazione } \end{cases} \]
Per farlo si utilizzerà sempre un test t, perchè non conoscendo la deviazione standard della popolazione, la curva t di Student userà quella del campione e incorporerà l’incertezza di questa approssimazione.
Per ottenere il dato medio di peso e lunghezza dei neonati, che sarà il nostro dato di popolazione, si è attinto al sito dell’ospedale pediatrico Bambino Gesù per cui risulta:
peso medio: 3300 grammi
lunghezza media: 50 cm
t.test(Peso, alternative = "two.sided", mu = 3300)
One Sample t-test
data: Peso
t = -1.505, df = 2497, p-value = 0.1324
alternative hypothesis: true mean is not equal to 3300
95 percent confidence interval:
3263.577 3304.791
sample estimates:
mean of x
3284.184 Il p-value supera di poco il 13% ed essendo quindi maggiore del livello di significatività, possiamo accettare l’ipotesi nulla: la media del campione è allineata con quella della popolazione
t.test(Lunghezza, alternative = "two.sided", mu = 500)
One Sample t-test
data: Lunghezza
t = -10.069, df = 2497, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 500
95 percent confidence interval:
493.6628 495.7287
sample estimates:
mean of x
494.6958 Lo stesso non può essere detto della lunghezza, per cui il p-value è molto vicino allo 0 e dobbiamo pertanto rifiutare l’ipotesi nulla: le lunghezze di questo campione non sono allineate con la media della popolazione.
3.3 Parti naturali vs cesarei
Concentriamoci innanzitutto sul descrivere queste variabili all’interno del nostro campione.
Tipo.parto
Ces Nat
728 1770 Il dataset mostra una netta maggioranza di parti naturali che diventa ancora più ovvia con l’ausilio di un grafico a barre
Provando ad ampliare la vista includendo ora anche gli ospedali abbiamo questo grafico:
Da questo grafico capiamo che nel campione sembra che gli ospedali seguano lo stesso andamento sui parti, quasi fossero equidistribuiti sui 3. Verrà saggiata questa ipotesi con dei test statistici ad hoc.
3.3.1 Si fanno più parti naturali che cesarei?
Il grafico sembrerebbe assolutamente indicare di si, ma possiamo essere così certi solo per quanto riguarda il nostro campione. Se volessimo generalizzare questa affermazione dovremmo verificare le probabilità di aver pescato un campione come il nostro.
Il test più adeguato per i confronti tra variabili categoriche è il Chi-quadrato. Può essere usato per verificare l’indipendenza di due variabili, o per verificare se una variabile si scosti significativamente da una distribuzione attesa (goodness of fit).
È proprio con questa seconda accezione che effettueremo il test, assumendo che non si facciano più parti di un tipo che di un altro. Quindi:
\[ \begin{cases} H_{0} : P\left( Cesareo \right) = 0,5 & \\ H_{1} : P\left( Cesareo \right) \neq 0,5\end{cases} \]
chisq.test(table(Tipo.parto))
Chi-squared test for given probabilities
data: table(Tipo.parto)
X-squared = 434.65, df = 1, p-value < 2.2e-16Il p-value è molto vicino allo 0, dobbiamo quindi rifiutare l’ipotesi nulla e possiamo affermare che effettivamente sussiste uno sbilanciamento nei tipi di parto, a favore dei naturali piuttosto che dei cesarei.
3.3.2 In alcuni ospedali si fanno più cesarei?
Dal momento che abbiamo dimostrato che non si fanno più cesarei che naturali, la domanda è da intendersi: ” C’è un ospedale che fa più cesarei degli altri? ”
Il grafico Figure 7 mostra che tra gli ospedali sembra esserci un’equidistribuzione dei cesarei, ma lasceremo al test saggiare questa ipotesi.
PartixOsp = neonati.clean %>%
dplyr::select(Ospedale, Tipo.parto)
PartixOsp %>%
filter(Tipo.parto == "Ces") %>%
mutate(Tipo.parto = droplevels(Tipo.parto)) %>%
table() %>%
chisq.test()
Chi-squared test for given probabilities
data: .
X-squared = 1, df = 2, p-value = 0.6065Il p-value è molto alto, poco oltre il 60%. Non rifiutiamo pertanto l’ipotesi nulla e affermiamo che tra gli ospedali c’è una distribuzione equa dei cesarei. Nessun ospedale ne fa quindi più degli altri in maniera statisticamente significativa.
3.3.3 Indipendenza delle variabili Ospedale e Tipo.parto
Infine, verifichiamo se le due variabili Ospedale e Tipo.parto siano indipendenti. La loro indipendenza implicherebbe che, conoscere il valore di una non consentirebbe di ricavare il valore dell’altra.
Il sistema di ipotesi pertanto è:
\[ \begin{cases} H_{0} : \mbox {Tipo.parto e Ospedale sono indipendenti} & \\ H_{1} : \mbox{Tipo.parto e Ospedale non sono indipendenti}\end{cases} \]
PartixOsp %>%
table() %>%
chisq.test()
Pearson's Chi-squared test
data: .
X-squared = 1.083, df = 2, p-value = 0.5819Il p-value è quasi al 60%, molto alto. Possiamo quindi accettare l’ipotesi nulla e affermare che le variabili Tipo.parto e Ospedale sono indipendenti. Questo vuol dire che non osserviamo nessuno schema ricorrente per cui, conosciuto l’ospedale, possiamo dire quale tipo di parto sarà praticato e viceversa.
4 Il modello di regressione
In questo capitolo verrà illustrato il procedimento seguito per definire un modello lineare basato sulle variabili del dataset e che consenta di prevedere in anticipo il peso del neonato.
Verrà seguita la procedura iterativa stepwise, partendo dal modello più complesso, con un regressore per ogni variabile del dataset e da lì si procederà per scremature sequenziali fino ad ottenere un buon compromesso tra varianza e semplicità.
I test e i criteri di valutazione che accompagneranno la procedura saranno:
- il t test sui coefficienti di regressione delle variabili e relativo p-value. Il sistema di ipotesi è il seguente
\[ \begin{cases} H_{0} : \beta_{i} = 0 &\\ H_{0} : \beta_{i} \neq 0 \end{cases} \]
Il coefficiente di determinazione aggiustato R^2 di ogni modello che indica la percentuale di adattabilità del modello ai dati veri
Il criterio di informazione Bayesiano (BIC)
Il test ANOVA
4.1 Verifiche preliminari
Il dataset si presenta con la variabile Fumatrici espressa con valori binari. Per facilitare l’analisi se ne creerà una nuova con valori “Si” e “No”.
detach(neonati.clean)
neonati.clean = neonati.clean %>%
mutate(Fumatrici_label = if_else(
Fumatrici==1, "Si", "No"
), Fumatrici_label = as.factor(Fumatrici_label)) %>%
dplyr::select(Anni.madre, N.gravidanze, Fumatrici, Fumatrici_label, everything())
attach(neonati.clean)
head(neonati.clean)| Anni.madre | N.gravidanze | Fumatrici | Fumatrici_label | Gestazione | Peso | Lunghezza | Cranio | Tipo.parto | Ospedale | Sesso |
|---|---|---|---|---|---|---|---|---|---|---|
| 26 | 0 | 0 | No | 42 | 3380 | 490 | 325 | Nat | osp3 | M |
| 21 | 2 | 0 | No | 39 | 3150 | 490 | 345 | Nat | osp1 | F |
| 34 | 3 | 0 | No | 38 | 3640 | 500 | 375 | Nat | osp2 | M |
| 28 | 1 | 0 | No | 41 | 3690 | 515 | 365 | Nat | osp2 | M |
| 20 | 0 | 0 | No | 38 | 3700 | 480 | 335 | Nat | osp3 | F |
| 32 | 0 | 0 | No | 40 | 3200 | 495 | 340 | Nat | osp2 | F |
Verifichiamo ora che la variabile risposta Peso, oggetto dell’analisi, sia distribuita secondo una gaussiana.
La variabile peso sembra essere distribuita come una gaussiana, tranne che per la coda sinistra, decisamente più allungata rispetto a quella destra.
D’altronde, ci si doveva aspettare una distribuzione del genere, avendo nel campione, delle nascite avvenute con gestazioni inferiori alla norma e di conseguenza neonati con peso inferiore (cfr. Figure 2).
Conduciamo il test di Shapiro-Wilk che verifica l’aderenza dei dati ad una normale per ulteriore scrupolo.
shapiro.test(Peso)
Shapiro-Wilk normality test
data: Peso
W = 0.97068, p-value < 2.2e-16Il p-value è sostanzialmente 0 e questo ci porta a rifiutare l’ipotesi nulla di normalità dei dati. Va ricordato che la variabile peso presentava comunque una forma lepticurtica e con un’asimmetria verso destra ( Table 2 ) e questo probabilmente incide sul test.
Si decide comunque di tenere a mente questo risultato ma spostare ogni valutazione una volta conseguito il modello definitivo.
4.2 Relazioni tra regressori e risposta
In questo paragrafo metteremo su grafico le correlazioni tra tutte le variabili del nostro dataset, con l’obiettivo di individuare a colpo d’occhio eventuali correlazioni che ci possano fornire una guida in sede di definizione del modello.
Lasciando per ora da parte le variabili categoriche come Fumatrici, Sesso, Tipo.parto per cui questo non è il grafico migliore, per tutte le quantitative possiamo ottenere alcune informazioni:
Anni.madre e N.gravidanze sembrano non avere singolarmente una correlazione con il peso
Gestazione, lunghezza e cranio ovviamente si. In particolare sembra che la gestazione abbia un andamento non lineare rispetto al peso, con la curva che tende ad appiattirsi per alti valori. Sembra ricalcare quasi una curva logaritmica.
Per quanto riguarda invece le variabili categoriche:
Su tutti i boxplot possiamo notare un tratto comune: la presenza di valori molto estremi, potenziali outliars che verranno valutati in sede di modellazione. Tendenzialmente sono verso il basso e questo deve essere dovuto al fatto che nel campione è presente un certo quantitativo di parti avvenuti con gestazioni inferiori alle 37 settimane, per la precisione:
[1] 162Quanto al sesso, abbiamo già verificato in questo paragrafo Section 3.1.1 che le differenze fisiche tra i sessi esistono e sono statisticamente rilevanti.
La variabile Fumatrici è molto sospetta perchè sembra indicare una cosa assolutamente controintuitiva, vale a dire la presenza di molti outliars sulle mamme NON fumatrici.
A un’analisi più attenta però, possiamo scartare qualsiasi considerazione fatta su questa variabile, perchè all’interno del campione abbiamo un numero di mamme fumatrici decisamente basso per poter trarre qualsiasi conclusione.
Fumatrici_label
No Si
0.96 0.04 Cerchiamo invece di capire se ci siano differenze statisticamente significative tra peso e tipo parto o ospdedale.
t.test(Peso~Tipo.parto, data = neonati.clean, alternative = "two.sided")
Welch Two Sample t-test
data: Peso by Tipo.parto
t = -0.13626, df = 1494.4, p-value = 0.8916
alternative hypothesis: true difference in means between group Ces and group Nat is not equal to 0
95 percent confidence interval:
-46.44246 40.40931
sample estimates:
mean in group Ces mean in group Nat
3282.047 3285.063 Il test t conferma, con un p-value a 89% che possiamo accettare l’ipotesi nulla per cui non esiste significativa differenza di peso rispetto al tipo di parto e in effetti, la cosa ha senso.
neonati.clean %>%
filter(Ospedale != "osp3") %>%
t.test(Peso~Ospedale, data= ., alternative = "two.sided")
Welch Two Sample t-test
data: Peso by Ospedale
t = -0.0092936, df = 1656.4, p-value = 0.9926
alternative hypothesis: true difference in means between group osp1 and group osp2 is not equal to 0
95 percent confidence interval:
-51.13390 50.65161
sample estimates:
mean in group osp1 mean in group osp2
3270.266 3270.507 neonati.clean %>%
filter(Ospedale != "osp2") %>%
t.test(Peso~Ospedale, data= ., alternative = "two.sided")
Welch Two Sample t-test
data: Peso by Ospedale
t = -1.6004, df = 1643.2, p-value = 0.1097
alternative hypothesis: true difference in means between group osp1 and group osp3 is not equal to 0
95 percent confidence interval:
-92.233560 9.348156
sample estimates:
mean in group osp1 mean in group osp3
3270.266 3311.709 neonati.clean %>%
filter(Ospedale != "osp1") %>%
t.test(Peso~Ospedale, data= ., alternative = "two.sided")
Welch Two Sample t-test
data: Peso by Ospedale
t = -1.623, df = 1680, p-value = 0.1048
alternative hypothesis: true difference in means between group osp2 and group osp3 is not equal to 0
95 percent confidence interval:
-90.993927 8.590812
sample estimates:
mean in group osp2 mean in group osp3
3270.507 3311.709 Anche nel confronto tra peso e ospedale il test t restituisce un p-value per tutte le 3 coppie (osp1-2, osp1-3, osp2-3) superiore al livello di significatività del 5%, ergo non rifiutiamo l’ipotesi nulla e concludiamo che la media dei pesi sugli ospedali non è diversa in maniera statisticamente significativa.
4.3 Definizione del modello con il metodo stepwise
Si partirà dal modello più complesso e si procederà per semplificazioni iterative, eliminando via via i regressori che non apportano un contenuto significativo.
4.3.1 mod 0
Partiamo dal modello che incorpora tutte le variabili, con l’unica eccezione della variabile Fumatrici che presenta la stessa informazione della variabile Fumatrici_label, con la differenza che la seconda è categorica mentre la prima è binaria.
mod0 = lm(Peso~.- Fumatrici, data = neonati.clean)
summary(mod0)
Call:
lm(formula = Peso ~ . - Fumatrici, data = neonati.clean)
Residuals:
Min 1Q Median 3Q Max
-1123.26 -181.53 -14.45 161.05 2611.89
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6735.7960 141.4790 -47.610 < 2e-16 ***
Anni.madre 0.8018 1.1467 0.699 0.4845
N.gravidanze 11.3812 4.6686 2.438 0.0148 *
Fumatrici_labelSi -30.2741 27.5492 -1.099 0.2719
Gestazione 32.5773 3.8208 8.526 < 2e-16 ***
Lunghezza 10.2922 0.3009 34.207 < 2e-16 ***
Cranio 10.4722 0.4263 24.567 < 2e-16 ***
Tipo.partoNat 29.6335 12.0905 2.451 0.0143 *
Ospedaleosp2 -11.0912 13.4471 -0.825 0.4096
Ospedaleosp3 28.2495 13.5054 2.092 0.0366 *
SessoM 77.5723 11.1865 6.934 5.18e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 274 on 2487 degrees of freedom
Multiple R-squared: 0.7289, Adjusted R-squared: 0.7278
F-statistic: 668.7 on 10 and 2487 DF, p-value: < 2.2e-16I dati salienti del modello che verranno utilizzati per tutti i confronti successivi sono:
il p- value dei coefficienti di regressione. In questo caso abbiamo che i coefficienti di regressioni significativamente vicini allo 0 sono:
gli anni della madre
le fumatrici
gli ospedali
R^2 adjusted pari a 0,728
Procediamo in maniera iterativa e proviamo eliminando le fumatrici, perchè effettivamente sono talmente poche rispetto al campione che non possono fornire rilievo nelle valutazioni.
4.3.2 mod 1
mod1 = update(mod0, ~.-Fumatrici_label)
summary(mod1)
Call:
lm(formula = Peso ~ Anni.madre + N.gravidanze + Gestazione +
Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso, data = neonati.clean)
Residuals:
Min 1Q Median 3Q Max
-1122.63 -181.96 -14.91 161.39 2615.07
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6735.4444 141.4845 -47.606 < 2e-16 ***
Anni.madre 0.8118 1.1467 0.708 0.4791
N.gravidanze 11.1201 4.6627 2.385 0.0172 *
Gestazione 32.3210 3.8138 8.475 < 2e-16 ***
Lunghezza 10.3064 0.3006 34.285 < 2e-16 ***
Cranio 10.4766 0.4263 24.577 < 2e-16 ***
Tipo.partoNat 29.3770 12.0888 2.430 0.0152 *
Ospedaleosp2 -11.0363 13.4475 -0.821 0.4119
Ospedaleosp3 28.5194 13.5038 2.112 0.0348 *
SessoM 77.3928 11.1858 6.919 5.77e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 274 on 2488 degrees of freedom
Multiple R-squared: 0.7288, Adjusted R-squared: 0.7278
F-statistic: 742.8 on 9 and 2488 DF, p-value: < 2.2e-16È rimasto tutto invariato, a conferma che effettivamente, la variabile Fumatrici non forniva un’informazione rilevante.
4.3.3 mod 2
La candidata successiva diventa la variabile Ospedale. Si sta aspettando ad eliminare quella che dai numeri sembrerebbe la prima in assoluto da scartare, gli anni della madre, per vedere se la situazione possa cambiare concentrandoci su variabili, come l’ospedale, concettualmente più separate dal peso del bambino.
mod2 = update(mod1, ~.-Ospedale)
summary(mod2)
Call:
lm(formula = Peso ~ Anni.madre + N.gravidanze + Gestazione +
Lunghezza + Cranio + Tipo.parto + Sesso, data = neonati.clean)
Residuals:
Min 1Q Median 3Q Max
-1139.50 -181.60 -14.59 160.14 2633.16
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6737.9269 141.5747 -47.593 < 2e-16 ***
Anni.madre 0.8793 1.1479 0.766 0.4438
N.gravidanze 11.4176 4.6676 2.446 0.0145 *
Gestazione 32.6300 3.8180 8.546 < 2e-16 ***
Lunghezza 10.2839 0.3009 34.176 < 2e-16 ***
Cranio 10.4896 0.4268 24.574 < 2e-16 ***
Tipo.partoNat 30.1222 12.1038 2.489 0.0129 *
SessoM 77.8374 11.2008 6.949 4.67e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 274.4 on 2490 degrees of freedom
Multiple R-squared: 0.7278, Adjusted R-squared: 0.727
F-statistic: 950.9 on 7 and 2490 DF, p-value: < 2.2e-16Abbiamo perso circa 1 millesimo su R^2 adjusted, ora a 0,727, assolutamente accettabile e conferma definitiva che la variabile Ospedale non apportasse significativo contributo al modello.
4.3.4 mod 3
A questo punto, non resta che provare a scartare gli anni della madre e vedere che impatto otteniamo. Ce lo si aspetta contenuto, essendo quella con il coefficiente di regressione con p-value maggiore.
mod3 = update(mod2, ~.-Anni.madre)
summary(mod3)
Call:
lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
Tipo.parto + Sesso, data = neonati.clean)
Residuals:
Min 1Q Median 3Q Max
-1129.14 -181.97 -16.26 160.95 2638.18
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6708.0171 136.0715 -49.298 < 2e-16 ***
N.gravidanze 12.7356 4.3385 2.935 0.00336 **
Gestazione 32.3253 3.7969 8.514 < 2e-16 ***
Lunghezza 10.2833 0.3009 34.177 < 2e-16 ***
Cranio 10.5063 0.4263 24.648 < 2e-16 ***
Tipo.partoNat 30.1601 12.1027 2.492 0.01277 *
SessoM 77.9171 11.1994 6.957 4.42e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 274.4 on 2491 degrees of freedom
Multiple R-squared: 0.7277, Adjusted R-squared: 0.727
F-statistic: 1109 on 6 and 2491 DF, p-value: < 2.2e-16I sospetti vengono confermati: R^2 adjusted rimane fermo a quota 0,727. Tutti i coefficienti di regressione sono ora sotto la soglia del 5% di livello di significatività, laddove il più alto p-value è a 1,2% per il tipo di parto.
4.3.5 mod 4
Proviamo ad eliminare la variabile tipo parto e vediamo cosa accade al modello. In questo caso ci aspetteremo di vedere R^2 cambiare.
mod4 = update(mod3, ~.-Tipo.parto)
summary(mod4)
Call:
lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
Sesso, data = neonati.clean)
Residuals:
Min 1Q Median 3Q Max
-1149.37 -180.98 -15.57 163.69 2639.09
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6681.7251 135.8036 -49.201 < 2e-16 ***
N.gravidanze 12.4554 4.3416 2.869 0.00415 **
Gestazione 32.3827 3.8008 8.520 < 2e-16 ***
Lunghezza 10.2455 0.3008 34.059 < 2e-16 ***
Cranio 10.5410 0.4265 24.717 < 2e-16 ***
SessoM 77.9807 11.2111 6.956 4.47e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 274.7 on 2492 degrees of freedom
Multiple R-squared: 0.727, Adjusted R-squared: 0.7265
F-statistic: 1327 on 5 and 2492 DF, p-value: < 2.2e-16Si perde un ulteriore millesimo su R^2 adjusted, ora a quota 0,726. Si ritiene la perdita accettabile in cambio di una semplificazione del modello. Ad ora infatti sono state scartate 4 variabili su 8 a fronte di una perdita di 2 millesimi su R^2 adjusted.
4.3.6 mod 5
Proviamo a scartare il numero di gravidanze pregresse per valutare cosa accade al modello. È una variabile significativa, ma è quella che, pur con p-value bassissimo, lo ha ben più alto di tutte le altre.
mod5 = update(mod4, ~.-N.gravidanze)
summary(mod5)
Call:
lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso,
data = neonati.clean)
Residuals:
Min 1Q Median 3Q Max
-1138.1 -184.4 -17.4 163.3 2626.7
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6651.6732 135.5952 -49.055 < 2e-16 ***
Gestazione 31.3262 3.7884 8.269 < 2e-16 ***
Lunghezza 10.2024 0.3009 33.909 < 2e-16 ***
Cranio 10.6706 0.4247 25.126 < 2e-16 ***
SessoM 79.1027 11.2205 7.050 2.31e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 275.1 on 2493 degrees of freedom
Multiple R-squared: 0.7261, Adjusted R-squared: 0.7257
F-statistic: 1652 on 4 and 2493 DF, p-value: < 2.2e-16R^2 adjusted rimane in effetti invariato e giustificherebbe l’eliminazione della variabile dal modello. Dal momento che gli anni della madre sono comunque ritenuti, in ambito medico, un fattore che incide sulle gravidanze, si effettueranno ulteriori test per capire se mantenere o meno la variabile.
4.3.7 Test sui modelli
Innanzitutto isoliamo le valutazioni agli ultimi 2, il modello 4 e il 5.
Il primo test che effettuiamo è l’ANOVA, acronimo di analysis of variances. Il test confronta il tasso di variazione degli scarti quadratici nel passaggio da un modello più complesso a uno più semplice, pesando questa variazione rispetto a quanto il modello più complesso riusciva a spiegare.
anova(mod5, mod4)| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 2493 | 188663107 | NA | NA | NA | NA |
| 2492 | 188042054 | 1 | 621053.3 | 8.230419 | 0.0041541 |
Il p-value viene molto basso, questo porterebbe a concludere che l’ipotesi nulla per cui il modello più semplice non fa perdere informazioni rispetto a quello più complesso debba essere rifiutata. L’ANOVA dunque suggerisce di mantenere il modello 4.
Proviamo con il test di informazione bayesiano.
BIC(mod1, mod2, mod3, mod4, mod5)| df | BIC | |
|---|---|---|
| mod1 | 11 | 35208.84 |
| mod2 | 9 | 35202.49 |
| mod3 | 8 | 35195.25 |
| mod4 | 7 | 35193.65 |
| mod5 | 6 | 35194.06 |
Il test passa sul modello che totalizza il punteggio più basso e anche in questo caso è il modello 4 a essere preferito.
Proviamo infine a lanciare la procedura stepwise automatizzata per verificare quali passi segua e se atterri sulle nostre stesse conclusioni
step.mod = stepAIC(mod0,
direction = "both",
k=log(length(Peso)))Start: AIC=28118.61
Peso ~ (Anni.madre + N.gravidanze + Fumatrici + Fumatrici_label +
Gestazione + Lunghezza + Cranio + Tipo.parto + Ospedale +
Sesso) - Fumatrici
Df Sum of Sq RSS AIC
- Anni.madre 1 36710 186779904 28111
- Fumatrici_label 1 90677 186833870 28112
- Ospedale 2 687555 187430749 28112
- N.gravidanze 1 446244 187189438 28117
- Tipo.parto 1 451073 187194266 28117
<none> 186743194 28119
- Sesso 1 3610705 190353899 28159
- Gestazione 1 5458852 192202046 28183
- Cranio 1 45318506 232061700 28654
- Lunghezza 1 87861708 274604902 29074
Step: AIC=28111.28
Peso ~ N.gravidanze + Fumatrici_label + Gestazione + Lunghezza +
Cranio + Tipo.parto + Ospedale + Sesso
Df Sum of Sq RSS AIC
- Fumatrici_label 1 91599 186871503 28105
- Ospedale 2 693914 187473818 28105
- Tipo.parto 1 452049 187231953 28110
<none> 186779904 28111
- N.gravidanze 1 631082 187410986 28112
+ Anni.madre 1 36710 186743194 28119
- Sesso 1 3617809 190397713 28151
- Gestazione 1 5424800 192204704 28175
- Cranio 1 45569477 232349381 28649
- Lunghezza 1 87852027 274631931 29066
Step: AIC=28104.68
Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
Ospedale + Sesso
Df Sum of Sq RSS AIC
- Ospedale 2 702925 187574428 28098
- Tipo.parto 1 444404 187315907 28103
<none> 186871503 28105
- N.gravidanze 1 608136 187479640 28105
+ Fumatrici_label 1 91599 186779904 28111
+ Anni.madre 1 37633 186833870 28112
- Sesso 1 3601860 190473363 28144
- Gestazione 1 5358199 192229702 28168
- Cranio 1 45613331 232484834 28642
- Lunghezza 1 88259386 275130889 29063
Step: AIC=28098.41
Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
Sesso
Df Sum of Sq RSS AIC
- Tipo.parto 1 467626 188042054 28097
<none> 187574428 28098
- N.gravidanze 1 648873 188223301 28099
+ Ospedale 2 702925 186871503 28105
+ Fumatrici_label 1 100610 187473818 28105
+ Anni.madre 1 44184 187530244 28106
- Sesso 1 3644818 191219246 28139
- Gestazione 1 5457887 193032315 28162
- Cranio 1 45747094 233321522 28636
- Lunghezza 1 87955701 275530129 29051
Step: AIC=28096.81
Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
Df Sum of Sq RSS AIC
<none> 188042054 28097
- N.gravidanze 1 621053 188663107 28097
+ Tipo.parto 1 467626 187574428 28098
+ Ospedale 2 726146 187315907 28103
+ Fumatrici_label 1 92548 187949505 28103
+ Anni.madre 1 45366 187996688 28104
- Sesso 1 3650790 191692844 28137
- Gestazione 1 5477493 193519547 28161
- Cranio 1 46098547 234140601 28637
- Lunghezza 1 87532691 275574744 29044summary(step.mod)
Call:
lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
Sesso, data = neonati.clean)
Residuals:
Min 1Q Median 3Q Max
-1149.37 -180.98 -15.57 163.69 2639.09
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6681.7251 135.8036 -49.201 < 2e-16 ***
N.gravidanze 12.4554 4.3416 2.869 0.00415 **
Gestazione 32.3827 3.8008 8.520 < 2e-16 ***
Lunghezza 10.2455 0.3008 34.059 < 2e-16 ***
Cranio 10.5410 0.4265 24.717 < 2e-16 ***
SessoM 77.9807 11.2111 6.956 4.47e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 274.7 on 2492 degrees of freedom
Multiple R-squared: 0.727, Adjusted R-squared: 0.7265
F-statistic: 1327 on 5 and 2492 DF, p-value: < 2.2e-16Anche la procedura automatizzata atterra sul modello 4 che a questo punto teniamo come quello migliore finora.
4.3.8 Relazioni non lineari dei regressori
Osservando l’immagine Figure 9, possiamo notare come la relazione tra la variabile risposta Peso e le variabili indipendenti Gestazione, Lunghezza e Cranio, non sia proprio lineare.
Si può notare infatti un cambio di pendenza nella retta di correlazione tra le variabili, che lascia un’apertura nei confronti di relazioni delle variabili indipendenti con la variabile risposta non lineari.
Se ne indagheranno alcune per verificare se si riesca a ottenere un miglioramento del R^2 adjusted rispetto al modello 4.
4.3.8.1 Risposta logaritmica della Gestazione
L’andamento del peso rispetto alla gestazione potrebbe far intendere una relazione logaritmica. Si prova quindi ad aggiornare il modello 4 con il logaritmo della gestazione.
mod6 = lm(formula = Peso ~ N.gravidanze + I(log(Gestazione)) + Lunghezza + Cranio +
Sesso, data = neonati.clean)
summary(mod6)
Call:
lm(formula = Peso ~ N.gravidanze + I(log(Gestazione)) + Lunghezza +
Cranio + Sesso, data = neonati.clean)
Residuals:
Min 1Q Median 3Q Max
-1153.33 -180.96 -15.52 162.10 2638.06
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -9680.4831 440.1171 -21.995 < 2e-16 ***
N.gravidanze 12.2975 4.3449 2.830 0.00469 **
I(log(Gestazione)) 1163.2891 140.9956 8.251 2.52e-16 ***
Lunghezza 10.2554 0.3026 33.886 < 2e-16 ***
Cranio 10.5297 0.4274 24.639 < 2e-16 ***
SessoM 78.7070 11.2197 7.015 2.95e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 274.9 on 2492 degrees of freedom
Multiple R-squared: 0.7265, Adjusted R-squared: 0.726
F-statistic: 1324 on 5 and 2492 DF, p-value: < 2.2e-16Nessun miglioramento rispetto al semplice modello 4, che viene confermato anche dal BIC.
BIC(mod4, mod6)| df | BIC | |
|---|---|---|
| mod4 | 7 | 35193.65 |
| mod6 | 7 | 35198.05 |
Si scarta quindi l’ipotesi di un andamento logaritmico.
4.3.8.2 Risposta quadratica rispetto alla lunghezza
Proviamo a sondare l’ipotesi che la variabile risposta sia meglio descritta con un andamento quadratico rispetto alla lunghezza.
mod7 = lm(formula = Peso ~ N.gravidanze + Gestazione + I(Lunghezza^2) + Cranio +
Sesso, data = neonati.clean)
summary(mod7)
Call:
lm(formula = Peso ~ N.gravidanze + Gestazione + I(Lunghezza^2) +
Cranio + Sesso, data = neonati.clean)
Residuals:
Min 1Q Median 3Q Max
-1161.44 -180.02 -11.17 165.90 2381.94
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -4.335e+03 1.442e+02 -30.059 < 2e-16 ***
N.gravidanze 1.314e+01 4.298e+00 3.057 0.00226 **
Gestazione 3.481e+01 3.713e+00 9.374 < 2e-16 ***
I(Lunghezza^2) 1.081e-02 3.075e-04 35.160 < 2e-16 ***
Cranio 1.047e+01 4.212e-01 24.845 < 2e-16 ***
SessoM 7.455e+01 1.110e+01 6.714 2.33e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 271.9 on 2492 degrees of freedom
Multiple R-squared: 0.7326, Adjusted R-squared: 0.7321
F-statistic: 1365 on 5 and 2492 DF, p-value: < 2.2e-16Abbiamo un incremento di ben un punto percentuale su R^2 adjusted che ora passa a 0,732, persino migliore del modello iniziale con tutte le variabili.
Verifichiamo il BIC rispetto al modello 4.
BIC(mod4, mod7)| df | BIC | |
|---|---|---|
| mod4 | 7 | 35193.65 |
| mod7 | 7 | 35142.06 |
Il modello 7 scalza il modello 4 e diventa il nostro modello migliore.
4.3.8.3 Risposta quadratica rispetto al Cranio
Come fatto per la lunghezza, testiamo una relazione quadratica rispetto al Cranio, partendo però non più dal modello 4, ma dal 7.
mod8 = lm(formula = Peso ~ N.gravidanze + Gestazione + I(Lunghezza^2) + I(Cranio^2) +
Sesso, data = neonati.clean)
summary(mod8)
Call:
lm(formula = Peso ~ N.gravidanze + Gestazione + I(Lunghezza^2) +
I(Cranio^2) + Sesso, data = neonati.clean)
Residuals:
Min 1Q Median 3Q Max
-1157.94 -179.42 -12.32 165.84 2366.68
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.653e+03 1.176e+02 -22.570 < 2e-16 ***
N.gravidanze 1.313e+01 4.289e+00 3.062 0.00222 **
Gestazione 3.651e+01 3.695e+00 9.882 < 2e-16 ***
I(Lunghezza^2) 1.083e-02 3.062e-04 35.360 < 2e-16 ***
I(Cranio^2) 1.560e-02 6.218e-04 25.081 < 2e-16 ***
SessoM 7.338e+01 1.108e+01 6.620 4.38e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 271.4 on 2492 degrees of freedom
Multiple R-squared: 0.7336, Adjusted R-squared: 0.7331
F-statistic: 1372 on 5 and 2492 DF, p-value: < 2.2e-16Un miglioramento di un ulteriore millesimo di punto su R^2 adjusted che porterebbe a ritenere che il modello 8 possa sostituire il 7.
Facciamo il test BIC per averne conferma.
BIC(mod4, mod7, mod8)| df | BIC | |
|---|---|---|
| mod4 | 7 | 35193.65 |
| mod7 | 7 | 35142.06 |
| mod8 | 7 | 35132.63 |
Abbiamo un miglioramento ulteriore rispetto al modello 7 confermato anche dal BIC e pertanto teniamo il modello 8 come migliore fino ad ora.
4.3.8.4 Risposta quadratica della Gestazione
Proviamo a saggiare la possibilità di una relazione non logaritmica, ma quadratica tra la Gustazione e la variabile risposta Peso.
mod9 = lm(formula = Peso ~ N.gravidanze + I(Gestazione^2) + I(Lunghezza^2) + I(Cranio^2) + Sesso,
data = neonati.clean)
summary(mod9)
Call:
lm(formula = Peso ~ N.gravidanze + I(Gestazione^2) + I(Lunghezza^2) +
I(Cranio^2) + Sesso, data = neonati.clean)
Residuals:
Min 1Q Median 3Q Max
-1156.66 -179.68 -12.29 166.35 2373.28
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.982e+03 7.042e+01 -28.145 < 2e-16 ***
N.gravidanze 1.313e+01 4.291e+00 3.059 0.00224 **
I(Gestazione^2) 4.840e-01 4.934e-02 9.810 < 2e-16 ***
I(Lunghezza^2) 1.087e-02 3.049e-04 35.635 < 2e-16 ***
I(Cranio^2) 1.565e-02 6.214e-04 25.179 < 2e-16 ***
SessoM 7.274e+01 1.109e+01 6.560 6.52e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 271.4 on 2492 degrees of freedom
Multiple R-squared: 0.7335, Adjusted R-squared: 0.7329
F-statistic: 1371 on 5 and 2492 DF, p-value: < 2.2e-16Vediamo che R^2 adjusted è rimasto invariato e secondo il principio del rasoio di Ockam, potremmo concludere che sia meglio mantenere il modello più semplice dal momento che non otteniamo risultati degni di nota con il quadrato della gestazione.
Effettuiamo un test bayesiano per avere ulteriore conferma
BIC(mod4, mod7, mod8, mod9)| df | BIC | |
|---|---|---|
| mod4 | 7 | 35193.65 |
| mod7 | 7 | 35142.06 |
| mod8 | 7 | 35132.63 |
| mod9 | 7 | 35134.00 |
Il test mostra che il modello 9, che coincide con il modello 8 con l’aggiunta della risposta quadratica sulla gestazione, non è da preferirsi al modello 8.
In linea con il test bayesiano e con il principio del rasoio di Ockam, si mantiene il modello 8 come il migliore ottenuto ad ora.
4.3.9 Multicolinearità
Ora che abbiamo identificato un modello definitivo, verifichiamo che non sussista multicolinearità tra le variabili indipendenti.
Se ci fossero infatti variabili correlate tra loro, potremmo avere effetti imprevedibili sulla risposta e quindi distorsioni nelle previsioni.
vif(mod = mod8) N.gravidanze Gestazione I(Lunghezza^2) I(Cranio^2) Sesso
1.023665 1.616818 1.996655 1.579492 1.041953 La soglia comunemente utilizzata è di 5. Dal momento che tutte le variabili sono abbondantemente sotto soglia, possiamo concludere che non sussiste multicolinearità tra loro.
4.3.10 Mean Squared Error (MSE)
Verifichiamo ora che l’errore (o residuo) medio quadratico sia il minimo possibile. Il confronto lo faremo ovviamente tra il nostro modello definitivo, il modello 4 e quello iniziale.
MSE_mod0 = mean(mod0$residuals^2)
MSE_mod4 = mean(mod4$residuals^2)
MSE_mod8 = mean(mod8$residuals^2)
MSE_mod0[1] 74757.08MSE_mod4[1] 75277.04MSE_mod8[1] 73460.61Decisamente il modello 8 è quello da selezionare. Pur essendo più semplice, vale a dire avendo meno variabili, cosa che porta tendenzialmente a un aumento dei residui, abbiamo addirittura una riduzione rispetto al modello iniziale.
4.3.11 Commenti sul modello 8
Il modello 8 presenta dunque il giusto compromesso tra semplicità e attinenza ai dati, con relazioni sia lineari che quadratiche che nel complesso hanno portato R^2 adjusted ad essere addirittura migliore che il valore inizialmente assunto con il modello avente tutte le variabili del campione.
Questi sono i coefficienti di regressione “beta” associati ad ogni regressore
mod8$coefficients (Intercept) N.gravidanze Gestazione I(Lunghezza^2) I(Cranio^2)
-2.653154e+03 1.313235e+01 3.651154e+01 1.082768e-02 1.559616e-02
SessoM
7.338031e+01 Al netto dell’intercetta, possiamo notare come ciascuno abbia un contributo positivo sulla variabile risposta Peso. In parole povere, ogni incremento nei regressori porta un aumento nel peso.
La variabile più impattante è quella del sesso, per cui i maschi incidono con 73 grammi in più rispetto alle femmine.
Segue poi la gestazione per cui, per ogni settimana di gestazione in più, il bambino ottiene 36,5 grammi. Ovviamente questo dato è da leggersi “cum grano salis” perchè il numero di settimane di gestazione ha un suo valore naturale e non può aumentare a dismisura.
Anche il numero di gravidanze pregresse sembra avere un effetto positivo sulla risposta del Peso, laddove per ogni gravidanza in più che la mamma ha avuto, si ha un incremento di 13 grammi nel peso del bambino.
I coefficienti di regressione più piccoli sono legati alla lunghezza e al diametro del cranio. Ciò è dovuto dal fatto che queste due variabili crescono con il quadrato, e questa crescita più rapida viene bilanciata da coefficienti di regressione più bassi. In questi casi, determinare l’incremento di peso, per incremento unitario della variabile, è più complesso.
Per ogni incremento unitario di lunghezza o diametro craniale, il bambino acquisisce rispettivamente questo quantitativo di grammi in più
\[ \beta_{i}\left( 1 + 2X^{n-1}_{i} \right) \]
Dove \(X^{n-1}_{i}\) è il valore della lunghezza o diametro craniale precedente all’aumento unitario.
4.3.12 Verifica della parte erratica
Come ultimo controllo di validità del nostro modello, si deve verificare che i residui aderiscano a delle caratteristiche che li tengano “sotto controllo”:
distribuzione normale: i residui devono essere distribuiti normalmente, vale a dire addensati attorno allo 0
omoschedasticità: la varianza deve essere costante
non correlazione: i residui non devono essere tra loro correlati. Altrimenti si avrebbe una tendenza all’aumento o diminuzione che porterebbe il modello a divergere (come per l’omoschedasticità del resto)
Effettueremo sia una verifica grafica, che una numerica.
4.3.12.1 Correlazione dei residui
Il grafico seguente mostra, per ogni risposta del modello (asse x), quanto scarto (asse y) ci sia rispetto al corrispettivo valore reale nel dataset.
Sostanzialmente stiamo vedendo, su grafico, come siano distribuiti i residui del nostro modello. Quello che si vuole vedere è una dispersione sostanzialmente orizzontale. Non vogliamo che il modello abbia una correlazione, positiva o negativa che sia, con i residui.
Il grafico sembra confermare che non ci sia correlazione dei residui
Per avere una conferma numerica a quello che il grafico ci mostra, effettuiamo il test di Durbin-Watson che ha questo sistema di ipotesi. Fissiamo come sempre il livello di significatività sempre al 5%.
\[ \begin{cases} H_{0} : \mbox{Autocorrelazione dei residui}=0 \\ H_{1} : \mbox{Autocorrelazione dei residui }>0 \end{cases} \]
dwtest(mod8)
Durbin-Watson test
data: mod8
DW = 1.9518, p-value = 0.1139
alternative hypothesis: true autocorrelation is greater than 0Con un p-value di poco superiore all’11% non rifiutiamo l’ipotesi nulla e confermiamo che i residui del modello non siano autocorrelati.
4.3.12.2 Distribuzione normale dei residui
Con il successivo grafico vogliamo verificare invece che i residui seguano una distribuzione normale. Per farlo utilizziamo un grafico q-q, che sta per quantile-quantile. Viene divisa in quantili la distribuzione normale e si posizionano sia sull’asse x che y. In questo modo, se i residui seguiranno questa distribuzione, dovranno posizionarsi sulla bisettrice del quadrante.
Ci si aspetta di vedere degli scostamenti, soprattutto alle code, poichè il dataset contiene bambini nati prematuramente e alcuni nati dopo le 40 settimane.
Come supposto inizialmente, i residui del modello sembrano addensarsi bene sulla bisettrice tranne che sulle code, in particolare sulla coda destra della normale (o sulla parte il alto a destra della bisettrice del grafico q-q).
Anche in questo caso effettuiamo un test numerico per vefiricare quanto i residui aderiscano a una normale standard.
shapiro.test(mod8$residuals)
Shapiro-Wilk normality test
data: mod8$residuals
W = 0.97856, p-value < 2.2e-16L’ipotesi nulla del test di Shapiro-Wilk è che il vettore di dati presentato sia distribuito normalmente. Il p-value è sostanzialmente 0, quindi dobbiamo rifiutare l’ipotesi nulla: i nostri residui non aderiscono bene a una distribuzione normale standard.
Per ora teniamo a mente questo risultato e proseguiamo con l’analisi senza invalidare il modello. Tireremo a valle tutte le conclusioni.
4.3.12.3 Omoschedasticità dei residui
Terzo criterio che i nostri residui devono rispettare è quello di omoschedasticità: la varianza dei residui non deve mostre dei pattern ben precisi, ma essere sparpagliata orizzontalmente.
Sull’asse delle x troviamo le risposte del modello, mentre su quello delle y troviamo le radici quadrate dei valori assoluti dei residui “studentizzati”. Questa operazione viene fatta al fine di addensare i residui quanto più possibile per rendere più facile l’individuazione degli outliars.
Il grafico lascia intendere che ci sia una distribuzione dei residui abbastanza sparpagliata, ma si nota come la retta di tendenza blu sia leggermente inclinata positivamente. Effettuiamo un test numerico per eliminare ogni dubbio.
Il test di Breusch-Pagan ha questo sistema di ipotesi
\[ \begin{cases} H_{0} : \mbox{I residui sono omoshedastici } \\ H_{1} : \mbox{I residui sono eteroschedastici } \end{cases} \]
bptest(mod8)
studentized Breusch-Pagan test
data: mod8
BP = 58.144, df = 5, p-value = 2.938e-11Il p-value viene bassissimo e quindi non possiamo accettare l’ipotesi nulla. Dobbiamo concludere che i residui non hanno una varianza costante, sono eteroschedastici.
Anche in questo caso, teniamo a mente il risultato e proseguiamo con l’analisi.
4.3.12.4 Leva e outliars
L’ultimo tassello da verificare è la presenza di outliars e valori di leva. Sono entrambi valori che indicano uno scostamento significativo del modello rispetto ai dati reali: nel caso degli outliars lo scostamento è sulla variabile risposta, per la leva è sui regressori.
Iniziamo con un grafico più semplice che rappresenti la distribuzione dei residui studentizzati con soglie a +2 e -2.
Il grafico sopra riportato mostra in effetti dei valori dei residui studentizzati oltre le soglie impostate. Questi sono punti di attenzione ma che non necessariamente si traducono in outliars perchè potrebbero non avere un effetto “dannoso” nei confronti del modello. A un’ispezione grafica infatti si riscontra come, a eccezione di pochi, siano tutti abbastanza vicini alle soglie. Verificheremo il tutto più avanti con un test specifico.
Il grafico seguente contiene molte informazioni, ma una volta chiarite dovrebbe risultare di facile lettura:
asse x: i valori di leva, che più sono alti e più indicano dei regressori lontani rispetto al resto
asse y: residui studentizzati
soglia leverage: linea tratteggiata blu che indica la soglia superata la quale identifichiamo i punti con leva da tenere sotto controllo. Non è detto che si traducano in outliars, ma potrebbero.
\[ \mbox{Soglia leverage}= \frac{2p}{n} \]
linea rossa tratteggiata: indica semplicemente lo 0 attorno al quale vogliamo si addensino i nostri residui
dimensione e colore dei punti: la dimensione dei punti aumenta con l’aumentare della distanza di Cook per ciascuno, come anche il colore che invece tenderà sempre più al rosso. La distanza di Cook è un parametro associato ai residui che, se superiore a 0,5 indica un punto da monitorare, se superiore a 1 indica un outliar.
In questo caso abbiamo valori di leva che si posizionano oltre la soglia, pochi comunque rispetto alla maggioranza, ma uno solo di loro, in alto a destra come potenziale outliar.
Effettuiamo comunque un test numerico per identificare altri outliars che possano sfuggire a un’analisi visiva del grafico
outlierTest(model = mod8) rstudent unadjusted p-value Bonferroni p
1549 9.028926 3.3838e-19 8.4529e-16
155 5.041006 4.9636e-07 1.2399e-03
1305 4.843831 1.3517e-06 3.3765e-03
1397 -4.287229 1.8780e-05 4.6912e-02Vengono identificati 4 valori outliars, valori che hanno un effetto particolarmente pronunciato nei confronti del modello. Ad ogni modo, trattandosi di 4 valori su un campione di 2500 osservazioni, qualunque effetto possano avere sarà sicuramente ammortizzato dal resto del campione.
5 Previsioni del modello
Proviamo ad effettuare, alcune previsioni di test sfruttando il modello.
predict(mod8, newdata = data.frame(N.gravidanze = 2,
Gestazione = 41,
Lunghezza = 500,
Cranio = 330,
Sesso = "F")) 1
3275.426 predict(mod8, newdata = data.frame(N.gravidanze = 0,
Gestazione = 40,
Lunghezza = 510,
Cranio = 330,
Sesso = "M")) 1
3395.389 predict(mod8, newdata = data.frame(N.gravidanze = 1,
Gestazione = 40,
Lunghezza = 500,
Cranio = 330,
Sesso = "M")) 1
3299.162 predict(mod8, newdata = data.frame(N.gravidanze = 1,
Gestazione = 40,
Lunghezza = 480,
Cranio = 330,
Sesso = "M")) 1
3086.94 I test sono stati condotti su nascite di cui si conosce il peso e, tranne nel primo caso, in cui la bambina è nata pesando 2,8 kg, quasi mezzo chilo di differenza rispetto al modello, negli altri casi la corrispondenza è molto vicina alla realtà.
Proviamo infine una previsione con dati parziali. Supponiamo ad esempio di voler prevedere il peso di una bambina la cui madre è alla terza gravidanza e partorirà alla 39ma settimana.
In questo caso, mancando i dati di lunghezza e cranio, utilizzeremo le medie del nostro campione.
predict(mod8, newdata = data.frame(N.gravidanze = 2,
Gestazione = 39,
Lunghezza = mean(neonati.clean$Lunghezza),
Cranio = mean(neonati.clean$Cranio),
Sesso = "F")) 1
3250.079 6 Conclusioni
Possiamo concludere che il modello si adatta bene sui valori “normali” e nascite non particolari, mentre sulle code, parti prematuri o tardivi, o anche pesi eccessivamente alti o bassi pur in gestazioni normali, non garantisce ottime performance.
Le evidenze riscontrate sull’eteroschedasticità e sulla distribuzione normale dei residui si consiglia di non farle incidere sulla valutazione del modello. Effettivamente, è il migliore identificabile, rimanendo nella linearità, con i dati a disposizione.
Si suggerisce pertanto di avviare una fase di “valutazione sul campo” del modello, usandolo nei vari reparti dei 3 ospedali per fare previsioni sul peso alla nascita. Il peso effettivo verrà poi confrontato con la previsione e verranno tratte nuove conclusioni sulle performance a valle.
Sarà inoltre possibile, in questo modo, fare una valutazione dei benefici apportati in termini economici tramite l’efficientamento delle risorse. Se i benefici apportati saranno comunque superiori dei costi sostenuti per i casi di previsione errata, sarà comunque accettabile mantenere un modello di cui si conoscono esattamente i limiti.